function [TEC] = TECVert(h1, h2, perigee, lat, lon, data, tol)
%
% Vertical TEC Numercial Integration
%
%DESCRIPTION:
%This function implements the numerical integration of the Vertical TEC as
%described in [2.5.8.1.1].
%This code is better explained in [2] and is used when:
% |lat2 - lat1| < 1e-5
% |lon2 - lon1| < 1e-5
%so when P1 and P2 are too close to eachother.
%
%PROTOTYPE:
% [TEC] = TECVert(g1, g2, data, tol)
%
%--------------------------------------------------------------------------
% INPUTS:
%   h1         [1x1]       1st Integration Height    [km]
%   h2         [1x1]       2nd Integration Height    [km]
%   lat        [1x1]       Latitude                  [deg] (of P1 or P2)
%   lon        [1x1]       Longitude                 [deg] (of P1 or P2)
%   data       [---]       Problem Data              [strc] (see NOTES)       
%   tol        [1x1]       Tolerance for Int. Accur. [-] (optional)
%--------------------------------------------------------------------------
% OUTPUTS:
%   TEC        [1x1]       Total Electron Content    [TECU]
%--------------------------------------------------------------------------
%
%NOTES:
% - For the tollerance (accuracy) "tol" the adviced value is 0.010.
% - The input "data" has been chosen to be a structure (for compactness of
%   the code) defined as:
%       data.a0      = 1st Az Coeff   [-]
%       data.a1      = 2nd Az Coeff   [-]
%       data.a2      = 3rd Az Coeff   [-]
%       data.mth     = Month          [month]
%       data.UT      = Universal Time [hours]
%       data.stModip = MODIP Table    [-] (modipneqg_wrapped.asc)
%       data.F2      = F2 Table       [-] (ccir21.asc)
%       data.Fm3     = Fm3 Table      [-] (ccir21.asc)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Constants
RE = 6371.2;                %Earth Mean Radius [km]

%Accuracy Input Check
if nargin == 3
    tol = 1e-2;
end

%Nr. of Discretization Points (Initial Guess)
n = 8;

%Integration
j = 0;          %First Run Index
res = 1;        %Residual for 1st run
while (res > tol)
    if j ~= 0
        %Doulbling Number of Points (from 2nd run on)
        n = 2*n;
    end
    %Integration Intervals
    Dn = (h2 - h1) / n;
    g = 0.5773502691896 * Dn;
    y = h1 + (Dn - g)/2;
    
    %Summatory Cycles
    aux = 0;
    for i = 1 : (n-1)
        %Determination of Altitude Case (1st Node)
        s1 = y + i*Dn;
        [N1] = ElecDens(s1, perigee, data, lat, lon);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Determination of Altitude Case (2nd Node)
        s2 = y + i*Dn + g;
        [N2] = ElecDens(s2, perigee, data, lat, lon);
        
        %Computation of the i-th Summatory Iteration
        aux = aux + (N1 + N2);
    end
    GN2 = Dn/2 * aux;
    if j == 0
        res = 1; %1t run
    else         %From 2nd run on
        res = abs(GN1 - GN2)/abs(GN1);
    end
    GN1 = GN2;
    %Define that 1st run is over
    j = 1;
end

%Total Electron Content
TEC = (GN2 + (GN2 - GN1)/15) * 1e-13;

end